home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / demo4.i < prev    next >
Text File  |  1995-07-26  |  2KB  |  93 lines
  1. /*
  2.    DEMO4.I
  3.    Airfoil demo
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1995.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. func demo4(mono=)
  11. /* DOCUMENT demo4
  12.          or demo4, mono=1
  13.      solves for the flow past a 2D airfoil using Kutta-Jakowski theory.
  14.      The colors represent static pressure (darker is lower pressure),
  15.      red lines are streamlines.
  16.      Solutions for various angles of attack are shown by animation.
  17.      With the mono=1 keyword, only the streamlines are shown.  (On a
  18.      monochrome terminal, the pressure makes the streamlines invisible.)
  19.  */
  20. {
  21.   require, "movie.i";
  22.   local nx, ny;
  23.   imax= 50;
  24.   movie, display;
  25.   display, imax;
  26. }
  27.  
  28. func display(i)
  29. {
  30.   attack= 28.*double(imax-i)/(imax-1);
  31.   solve, attack, 2.0, 0.2;
  32.   if (i==1) {
  33.     extern cn, cx;
  34.     cn= 0.8*min(pressure);
  35.     cx= max(pressure);
  36.   }
  37.   limits, -2.8, 2.4, -2.4, 2.8;
  38.   plmesh, z.im,z.re;
  39.   if (!mono) plf, zncen(pressure), cmin=cn,cmax=cx;
  40.   plc, potential.im, marks=0,color="red", levs=span(-2.5,3.5,33);
  41.   plg,z.im(,1),z.re(,1), marks=0;
  42.   return i<imax;
  43. }
  44.  
  45. func solve(attack, chord, thick)
  46. {
  47.   extern w, z, potential, velocity, pressure; /* outputs */
  48.   attack*= pi/180.;
  49.   a= 0.5*chord;
  50.   w= get_mesh(a, attack);
  51.   /* the log(w) term adds just enough circulation to move the rear
  52.      stagnation point to the trailing edge -- the Kutta condition */
  53.   potential= jakowski(w, a) + 2i*a*sin(attack)*log(w);
  54.   dpdw= 1.-(a/w)^2 + 2i*a*sin(attack)/w;
  55.   emith= exp(-1i*attack);
  56.   a= a*emith;
  57.   b= thick*emith;  /* could add camber here too someday? */
  58.   w-= b;
  59.   a-= b;
  60.   z= jakowski(w, a);
  61.   dzdw= 1.-(a/w)^2;
  62.   velocity= conj(dpdw/dzdw);
  63.   pressure= 0.5*(1.0-abs(velocity)^2);
  64. }
  65.  
  66. func get_mesh(a, attack)
  67. {
  68.   /* get a mesh in the w-plane
  69.      -- the plane in which the airfoil is a circle
  70.      the mesh splits at the point which will become the trailing edge */
  71.   if (is_void(nx)) nx= 120;
  72.   if (is_void(ny)) ny= 30;
  73.   a= abs(a);
  74.   theta= span(1.e-9, 2*pi-1.e-9, nx) - attack;
  75.   r= a/span(1.0+1.e-9, 0.25, ny)(-,);
  76.   return r*exp(1i*theta);
  77. }
  78.  
  79. func jakowski(z, a)
  80. {
  81.   /* Jakowski transform - circle of radius a into slot of length 4a */
  82.   return z + a*a/z;
  83. }
  84.  
  85. func ijakowski(z, a)
  86. {
  87.   /* Inverse Jakowski - slot back into circle (unused here) */
  88.   z*= 0.5;
  89.   a= complex(a);
  90.   sgn= sign((z*conj(a)).re);
  91.   return z + sgn(z);
  92. }
  93.